Load in packages…

Add in our data and make “Year” a catagorical predictor

#Read in Raw Data
TurtleData <- read.csv("Data/TidyNests.csv")
TurtleData$Year <- as.factor(TurtleData$Year)
TurtleData$X <- NULL

Make some Histograms for Nest Abundances and False Crawl Abundances

#Visualize our Counts Data. We can see that its not normally distributed, thats okay because we knew it wouldnt be
#Count Data falls under the Poisson distribution, so this visual makes sense to see
hist(TurtleData$NestAbundance)

### We'll Use This Information Later

Test our predictor variables for correlation among scales

Mean.Dep <- TurtleData[, c("MEAN.Depth.250", "MEAN.Depth.500", 
                          "MEAN.Depth.750", "MEAN.Depth.1000")]
ggpairs(Mean.Dep)

SD.Dep <- TurtleData[, c("STD.Depth.250", "STD.Depth.500", 
                        "STD.Depth.750", "STD.Depth.1000")]
ggpairs(SD.Dep)

Mean.Slp <- TurtleData[, c("MEAN.Slope.ON", "MEAN.Slope.250", "MEAN.Slope.500", 
                          "MEAN.Slope.750", "MEAN.Slope.1000")]
ggpairs(Mean.Slp)

SD.Slp <- TurtleData[, c("STD.Slope.ON", "STD.Slope.250", "STD.Slope.500", 
                        "STD.Slope.750", "STD.Slope.1000")]
ggpairs(SD.Slp)

Dists <- TurtleData[, c("Distance.to.0m", "Distance.to.10m", 
                       "Distance.to.20m", "Distance.to.30m")]

ggpairs(Dists)

Unsuprisingly, we see corrilation among scales here for most of our data. This makes sense thinking about what we know about the ocean, and we do not want to use the same variable more than once anyway. So before we actually build our global model, we want to use AICc to see what scales these predictor variables are acting on our response variables. Lets do some scale selection! First lets build a bunch of mono-variable models.

options(na.action = "na.fail")

for (i in (2:25)) {TurtleData[,i] <- c(scale(TurtleData[,i], center = TRUE, scale = TRUE))}

NullModel <- glmer(NestAbundance ~ (1|Year), data = TurtleData, family = poisson)

DepthMean250 <- glmer(NestAbundance ~ (1|Year) + MEAN.Depth.250, data = TurtleData, family = poisson)

DepthSTD250 <- glmer(NestAbundance ~ (1|Year) + STD.Depth.250, data = TurtleData, family = poisson)

SlopeMean250 <- glmer(NestAbundance ~ (1|Year) + MEAN.Slope.250, data = TurtleData, family = poisson)

SlopeSTD250 <- glmer(NestAbundance ~ (1|Year) + STD.Slope.250, data = TurtleData, family = poisson)


DepthMean500 <- glmer(NestAbundance ~ (1|Year) + MEAN.Depth.500, data = TurtleData, family = poisson)

DepthSTD500 <- glmer(NestAbundance ~ (1|Year) + STD.Depth.500, data = TurtleData, family = poisson)

SlopeMean500 <- glmer(NestAbundance ~ (1|Year) + MEAN.Slope.500, data = TurtleData, family = poisson)

SlopeSTD500 <- glmer(NestAbundance ~ (1|Year) + STD.Slope.500, data = TurtleData, family = poisson)


DepthMean750 <- glmer(NestAbundance ~ (1|Year) + MEAN.Depth.750, data = TurtleData, family = poisson)

DepthSTD750 <- glmer(NestAbundance ~ (1|Year) + STD.Depth.750, data = TurtleData, family = poisson)

SlopeMean750 <- glmer(NestAbundance ~ (1|Year) + MEAN.Slope.750, data = TurtleData, family = poisson)

SlopeSTD750 <- glmer(NestAbundance ~ (1|Year) + STD.Slope.750, data = TurtleData, family = poisson)


Rugosity250 <- glmer(NestAbundance ~ (1|Year) + STD.Depth.250 * STD.Slope.250, data = TurtleData, family = poisson)

Rugosity500 <- glmer(NestAbundance ~ (1|Year) + STD.Depth.500 * STD.Slope.500, data = TurtleData, family = poisson)

Rugosity750 <- glmer(NestAbundance ~ (1|Year) + STD.Depth.750 * STD.Slope.750, data = TurtleData, family = poisson)

Rugosity1000 <-glmer(NestAbundance ~ (1|Year) + STD.Depth.1000 * STD.Slope.1000, data = TurtleData, family = poisson)


DepthMean1000 <- glmer(NestAbundance ~ (1|Year) + MEAN.Depth.1000, data = TurtleData, family = poisson)

DepthSTD1000 <- glmer(NestAbundance ~ (1|Year) + STD.Depth.1000, data = TurtleData, family = poisson)

SlopeMean1000 <- glmer(NestAbundance ~ (1|Year) + MEAN.Slope.1000, data = TurtleData, family = poisson)

SlopeSTD1000 <- glmer(NestAbundance ~ (1|Year) + STD.Slope.1000, data = TurtleData, family = poisson)

Dist.to.0 <- glmer(NestAbundance ~ (1|Year) + Distance.to.0m, data = TurtleData, family = poisson)

Dist.to.10 <- glmer(NestAbundance ~ (1|Year) + Distance.to.10m, data = TurtleData, family = poisson)

Dist.to.20 <- glmer(NestAbundance ~ (1|Year) + Distance.to.20m, data = TurtleData, family = poisson)

Dist.to.30 <- glmer(NestAbundance ~ (1|Year) + Distance.to.30m, data = TurtleData, family = poisson)

Then, lets pit all the scale models at each variable against one another.

depth.mean.aic <- list("DepthMean250" = DepthMean250,  "DepthMean500" = DepthMean500,  
                       "DepthMean750" = DepthMean750, "DepthMean1000" = DepthMean1000, "Null" = NullModel)

model.sel(depth.mean.aic)
## Model selection table 
##               (Int) MEA.Dpt.250 MEA.Dpt.500 MEA.Dpt.750 MEA.Dpt.1000
## DepthMean1000 1.463                                          -0.2084
## DepthMean750  1.472                             -0.1631             
## DepthMean500  1.479                  -0.107                         
## Null          1.485                                                 
## DepthMean250  1.485   -0.006744                                     
##                     family df    logLik   AICc  delta weight
## DepthMean1000 poisson(log)  3 -2914.479 5835.0   0.00      1
## DepthMean750  poisson(log)  3 -2955.517 5917.1  82.08      0
## DepthMean500  poisson(log)  3 -2993.886 5993.8 158.81      0
## Null          poisson(log)  2 -3023.143 6050.3 215.32      0
## DepthMean250  poisson(log)  3 -3023.028 6052.1 217.10      0
## Models ranked by AICc(x) 
## Random terms (all models): 
## '1 | Year'
summary(DepthMean1000)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: NestAbundance ~ (1 | Year) + MEAN.Depth.1000
##    Data: TurtleData
## 
##      AIC      BIC   logLik deviance df.resid 
##   5835.0   5849.9  -2914.5   5829.0     1077 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7972 -1.1773 -0.2841  0.7750  7.0389 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Year   (Intercept) 0.1186   0.3444  
## Number of obs: 1080, groups:  Year, 10
## 
## Fixed effects:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      1.46336    0.10995   13.31   <2e-16 ***
## MEAN.Depth.1000 -0.20845    0.01421  -14.67   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## MEAN.D.1000 0.027
car::Anova(DepthMean1000)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: NestAbundance
##                 Chisq Df Pr(>Chisq)    
## MEAN.Depth.1000 215.1  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############

slope.mean.aic <- list("SlopeMean250" = SlopeMean250, "SlopeMean500" = SlopeMean500, 
                       "SlopeMean750" = SlopeMean750, "SlopeMean1000" = SlopeMean1000, "Null" = NullModel)

model.sel(slope.mean.aic)
## Model selection table 
##               (Int) MEA.Slp.250 MEA.Slp.500 MEA.Slp.750 MEA.Slp.1000
## SlopeMean250  1.457     -0.2367                                     
## SlopeMean500  1.465                 -0.2045                         
## SlopeMean750  1.463                             -0.2291             
## SlopeMean1000 1.466                                           -0.221
## Null          1.485                                                 
##                     family df    logLik   AICc  delta weight
## SlopeMean250  poisson(log)  3 -2885.812 5777.6   0.00      1
## SlopeMean500  poisson(log)  3 -2927.465 5861.0  83.31      0
## SlopeMean750  poisson(log)  3 -2929.919 5865.9  88.22      0
## SlopeMean1000 poisson(log)  3 -2943.952 5893.9 116.28      0
## Null          poisson(log)  2 -3023.143 6050.3 272.65      0
## Models ranked by AICc(x) 
## Random terms (all models): 
## '1 | Year'
summary(SlopeMean250)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: NestAbundance ~ (1 | Year) + MEAN.Slope.250
##    Data: TurtleData
## 
##      AIC      BIC   logLik deviance df.resid 
##   5777.6   5792.6  -2885.8   5771.6     1077 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1863 -1.1104 -0.2273  0.8743  5.6166 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Year   (Intercept) 0.1186   0.3444  
## Number of obs: 1080, groups:  Year, 10
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     1.45739    0.10997   13.25   <2e-16 ***
## MEAN.Slope.250 -0.23669    0.01442  -16.42   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## MEAN.Sl.250 0.030
car::Anova(SlopeMean250)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: NestAbundance
##                Chisq Df Pr(>Chisq)    
## MEAN.Slope.250 269.5  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############

rugosity.aic <- list("Rugosity250" = Rugosity250, "Rugosity500" = Rugosity500, 
                      "Rugosity750" = Rugosity750, "Rugosity1000" = Rugosity1000, "Null" = NullModel)

model.sel(rugosity.aic)
## Model selection table 
##              (Int) STD.Dpt.250 STD.Slp.250 STD.Dpt.250:STD.Slp.250 STD.Dpt.500
## Rugosity500  1.542                                                     0.01304
## Rugosity1000 1.502                                                            
## Rugosity750  1.485                                                            
## Rugosity250  1.398      0.0298     -0.2719                 0.07936            
## Null         1.485                                                            
##              STD.Slp.500 STD.Dpt.500:STD.Slp.500 STD.Dpt.750 STD.Slp.750
## Rugosity500      -0.2425                 -0.1725                        
## Rugosity1000                                                            
## Rugosity750                                          -0.2771       0.265
## Rugosity250                                                             
## Null                                                                    
##              STD.Dpt.750:STD.Slp.750 STD.Dpt.1000 STD.Slp.1000
## Rugosity500                                                   
## Rugosity1000                              -0.3159       0.3908
## Rugosity750                 -0.06816                          
## Rugosity250                                                   
## Null                                                          
##              STD.Dpt.1000:STD.Slp.1000       family df    logLik   AICc  delta
## Rugosity500                            poisson(log)  5 -2801.814 5613.7   0.00
## Rugosity1000                   -0.0898 poisson(log)  5 -2886.044 5782.1 168.46
## Rugosity750                            poisson(log)  5 -2888.796 5787.6 173.96
## Rugosity250                            poisson(log)  5 -2917.066 5844.2 230.50
## Null                                   poisson(log)  2 -3023.143 6050.3 436.61
##              weight
## Rugosity500       1
## Rugosity1000      0
## Rugosity750       0
## Rugosity250       0
## Null              0
## Models ranked by AICc(x) 
## Random terms (all models): 
## '1 | Year'
summary(Rugosity500)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: NestAbundance ~ (1 | Year) + STD.Depth.500 * STD.Slope.500
##    Data: TurtleData
## 
##      AIC      BIC   logLik deviance df.resid 
##   5613.6   5638.6  -2801.8   5603.6     1075 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0761 -1.0635 -0.2950  0.7714  6.0333 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Year   (Intercept) 0.1186   0.3444  
## Number of obs: 1080, groups:  Year, 10
## 
## Fixed effects:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  1.54203    0.11047  13.959  < 2e-16 ***
## STD.Depth.500                0.01304    0.02375   0.549    0.583    
## STD.Slope.500               -0.24246    0.03239  -7.486 7.08e-14 ***
## STD.Depth.500:STD.Slope.500 -0.17249    0.01903  -9.064  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) STD.Dp.500 STD.S.
## STD.Dpt.500 -0.014                  
## STD.Slp.500  0.068 -0.663           
## STD.D.500:S -0.088  0.116     -0.384
car::Anova(Rugosity500)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: NestAbundance
##                               Chisq Df Pr(>Chisq)    
## STD.Depth.500                 2.604  1     0.1066    
## STD.Slope.500               141.213  1     <2e-16 ***
## STD.Depth.500:STD.Slope.500  82.154  1     <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############

dist.to.dep.aic <- list("Dist.to.0" = Dist.to.0,
                        "Dist.to.10" = Dist.to.10,
                        "Dist.to.20" = Dist.to.20,
                        "Dist.to.30" = Dist.to.30)

model.sel(dist.to.dep.aic)
## Model selection table 
##            (Int) Dst.to.0m Dst.to.10m Dst.to.20m Dst.to.30m       family df
## Dist.to.20 1.429                          0.3331            poisson(log)  3
## Dist.to.30 1.433                                     0.3297 poisson(log)  3
## Dist.to.0  1.439   -0.3058                                  poisson(log)  3
## Dist.to.10 1.485              0.01826                       poisson(log)  3
##               logLik   AICc  delta weight
## Dist.to.20 -2737.011 5480.0   0.00      1
## Dist.to.30 -2774.551 5555.1  75.08      0
## Dist.to.0  -2793.227 5592.5 112.43      0
## Dist.to.10 -3022.285 6050.6 570.55      0
## Models ranked by AICc(x) 
## Random terms (all models): 
## '1 | Year'
summary(Dist.to.20)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: NestAbundance ~ (1 | Year) + Distance.to.20m
##    Data: TurtleData
## 
##      AIC      BIC   logLik deviance df.resid 
##     5480     5495    -2737     5474     1077 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9521 -1.0876 -0.2367  0.7884  5.3526 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Year   (Intercept) 0.1186   0.3444  
## Number of obs: 1080, groups:  Year, 10
## 
## Fixed effects:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      1.42862    0.11002   12.98   <2e-16 ***
## Distance.to.20m  0.33306    0.01398   23.83   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Dstnc.t.20m -0.043
car::Anova(Dist.to.20)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: NestAbundance
##                  Chisq Df Pr(>Chisq)    
## Distance.to.20m 567.73  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Now that we have selected the most relevant scale for each variable, build a new data frame and global model out of those variables at the most relevant scales to our response.

Turt.nest.data1 <- data.frame(nests = TurtleData$NestAbundance,
                         Year = as.numeric(TurtleData$Year),
                         muDep = abs(TurtleData$MEAN.Depth.1000),
                         muSlop = TurtleData$MEAN.Slope.250,
                         Rugosity = TurtleData$STD.Slope.500,
                         dist20 = TurtleData$Distance.to.20m,
                         shelfW = (TurtleData$Distance.to.20m - TurtleData$Distance.to.10m)
                         )

ggpairs(Turt.nest.data1, columns = 3:7) + theme_bw()

GlobalModel <- glmer(nests ~ muDep + muSlop + Rugosity + 
             dist20 + shelfW + (1|Year), data = Turt.nest.data1, family = poisson)

performance::check_collinearity(GlobalModel)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##      Term  VIF Increased SE Tolerance
##     muDep 1.15         1.07      0.87
##  Rugosity 1.33         1.16      0.75
##    shelfW 3.23         1.80      0.31
## 
## Moderate Correlation
## 
##    Term  VIF Increased SE Tolerance
##  muSlop 5.14         2.27      0.19
##  dist20 7.40         2.72      0.14
GlobalModel_Dist <- glmer(nests ~ muDep + muSlop + Rugosity + 
             dist20 + (1|Year), data = Turt.nest.data1, family = poisson)

performance::check_collinearity(GlobalModel_Dist)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##      Term  VIF Increased SE Tolerance
##     muDep 1.03         1.02      0.97
##    muSlop 2.41         1.55      0.41
##  Rugosity 1.30         1.14      0.77
##    dist20 2.79         1.67      0.36
GlobalModel_Shelf <- glmer(nests ~ muDep + muSlop + Rugosity + 
             shelfW + (1|Year), data = Turt.nest.data1, family = poisson)

performance::check_collinearity(GlobalModel_Shelf)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##      Term  VIF Increased SE Tolerance
##     muDep 1.08         1.04      0.93
##    muSlop 1.08         1.04      0.92
##  Rugosity 1.32         1.15      0.76
##    shelfW 1.22         1.11      0.82

For now for simplicities sake, we can dredge our global model (while ignoring our multicoliniar variables)

dredge(GlobalModel_Shelf)
## Fixed term is "(Intercept)"
## Global model call: glmer(formula = nests ~ muDep + muSlop + Rugosity + shelfW + 
##     (1 | Year), data = Turt.nest.data1, family = poisson)
## ---
## Model selection table 
##    (Intrc)   muDep   muSlp   Rgsty  shlfW df    logLik   AICc  delta weight
## 16   1.550 -0.1726 -0.1948 -0.1862 0.1990  6 -2638.878 5289.8   0.00      1
## 15   1.406         -0.2050 -0.1560 0.2149  5 -2654.294 5318.6  28.81      0
## 12   1.525 -0.1278 -0.2270         0.2547  5 -2670.354 5350.8  60.93      0
## 11   1.416         -0.2312         0.2612  4 -2679.250 5366.5  76.70      0
## 14   1.614 -0.2325         -0.2644 0.1869  5 -2717.685 5445.4 155.59      0
## 8    1.658 -0.2773 -0.1642 -0.3159         5 -2731.240 5472.5 182.70      0
## 13   1.421                 -0.2248 0.2087  4 -2745.887 5499.8 209.98      0
## 7    1.429         -0.1762 -0.2786         4 -2773.047 5554.1 264.30      0
## 10   1.589 -0.1780                 0.2782  4 -2778.305 5564.6 274.81      0
## 6    1.703 -0.3204         -0.3868         4 -2791.655 5591.3 301.51      0
## 9    1.438                         0.2866  3 -2795.565 5597.2 307.32      0
## 5    1.440                 -0.3463         3 -2847.544 5701.1 411.28      0
## 4    1.662 -0.2414 -0.2320                 4 -2852.347 5712.7 422.90      0
## 3    1.457         -0.2367                 3 -2885.812 5777.6 487.81      0
## 2    1.729 -0.2880                         3 -2976.166 5958.4 668.52      0
## 1    1.485                                 2 -3023.143 6050.3 760.46      0
## Models ranked by AICc(x) 
## Random terms (all models): 
## '1 | Year'
shelf.res <- dredge(GlobalModel_Shelf)
## Fixed term is "(Intercept)"
importance(shelf.res)
##                      muSlop shelfW Rugosity muDep
## Sum of weights:      1      1      1        1    
## N containing models: 8      8      8        8
car::Anova(GlobalModel_Shelf)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: nests
##            Chisq Df Pr(>Chisq)    
## muDep     30.729  1  2.967e-08 ***
## muSlop   153.033  1  < 2.2e-16 ***
## Rugosity  56.643  1  5.227e-14 ***
## shelfW   176.732  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

And finally, lets plot our pairwise variables for our best model

ggplot(Turt.nest.data1, aes(muDep, nests)) +
  geom_point(size = 1) +
   geom_smooth(method = "glm", se = FALSE,
       method.args = list(family = "poisson"))
## `geom_smooth()` using formula 'y ~ x'

ggplot(Turt.nest.data1, aes(muSlop, nests)) +
  geom_point(size = 1) +
   geom_smooth(method = "glm", se = FALSE,
       method.args = list(family = "poisson"))
## `geom_smooth()` using formula 'y ~ x'

ggplot(Turt.nest.data1, aes(Rugosity, nests)) +
  geom_point(size = 1) +
   geom_smooth(method = "glm", se = FALSE,
       method.args = list(family = "poisson"))
## `geom_smooth()` using formula 'y ~ x'

ggplot(Turt.nest.data1, aes(shelfW, nests)) +
  geom_point(size = 1) +
   geom_smooth(method = "glm", se = FALSE,
       method.args = list(family = "poisson"))
## `geom_smooth()` using formula 'y ~ x'

What does this look like biologically?

LETS DO IT ALL AGAIN, this time with False Crawls

#Read in Raw Data
TurtleData <- read.csv("Data/TidyFC.csv")
TurtleData$Year <- as.factor(TurtleData$Year)
TurtleData$X <- NULL
for (i in (2:25)) {TurtleData[,i] <- c(scale(TurtleData[,i], center = TRUE, scale = TRUE))}

NullModel <- glmer(FCAbundance ~ (1|Year), data = TurtleData, family = poisson)

DepthMean250 <- glmer(FCAbundance ~ (1|Year) + MEAN.Depth.250, data = TurtleData, family = poisson)

DepthSTD250 <- glmer(FCAbundance ~ (1|Year) + STD.Depth.250, data = TurtleData, family = poisson)

SlopeMean250 <- glmer(FCAbundance ~ (1|Year) + MEAN.Slope.250, data = TurtleData, family = poisson)

SlopeSTD250 <- glmer(FCAbundance ~ (1|Year) + STD.Slope.250, data = TurtleData, family = poisson)


DepthMean500 <- glmer(FCAbundance ~ (1|Year) + MEAN.Depth.500, data = TurtleData, family = poisson)

DepthSTD500 <- glmer(FCAbundance ~ (1|Year) + STD.Depth.500, data = TurtleData, family = poisson)

SlopeMean500 <- glmer(FCAbundance ~ (1|Year) + MEAN.Slope.500, data = TurtleData, family = poisson)

SlopeSTD500 <- glmer(FCAbundance ~ (1|Year) + STD.Slope.500, data = TurtleData, family = poisson)


DepthMean750 <- glmer(FCAbundance ~ (1|Year) + MEAN.Depth.750, data = TurtleData, family = poisson)

DepthSTD750 <- glmer(FCAbundance ~ (1|Year) + STD.Depth.750, data = TurtleData, family = poisson)

SlopeMean750 <- glmer(FCAbundance ~ (1|Year) + MEAN.Slope.750, data = TurtleData, family = poisson)

SlopeSTD750 <- glmer(FCAbundance ~ (1|Year) + STD.Slope.750, data = TurtleData, family = poisson)


Rugosity250 <- glmer(FCAbundance ~ (1|Year) + STD.Depth.250 * STD.Slope.250, data = TurtleData, family = poisson)

Rugosity500 <- glmer(FCAbundance ~ (1|Year) + STD.Depth.500 * STD.Slope.500, data = TurtleData, family = poisson)

Rugosity750 <- glmer(FCAbundance ~ (1|Year) + STD.Depth.750 * STD.Slope.750, data = TurtleData, family = poisson)

Rugosity1000 <-glmer(FCAbundance ~ (1|Year) + STD.Depth.1000 * STD.Slope.1000, data = TurtleData, family = poisson)


DepthMean1000 <- glmer(FCAbundance ~ (1|Year) + MEAN.Depth.1000, data = TurtleData, family = poisson)

DepthSTD1000 <- glmer(FCAbundance ~ (1|Year) + STD.Depth.1000, data = TurtleData, family = poisson)

SlopeMean1000 <- glmer(FCAbundance ~ (1|Year) + MEAN.Slope.1000, data = TurtleData, family = poisson)

SlopeSTD1000 <- glmer(FCAbundance ~ (1|Year) + STD.Slope.1000, data = TurtleData, family = poisson)

Dist.to.0 <- glmer(FCAbundance ~ (1|Year) + Distance.to.0m, data = TurtleData, family = poisson)

Dist.to.10 <- glmer(FCAbundance ~ (1|Year) + Distance.to.10m, data = TurtleData, family = poisson)

Dist.to.20 <- glmer(FCAbundance ~ (1|Year) + Distance.to.20m, data = TurtleData, family = poisson)

Dist.to.30 <- glmer(FCAbundance ~ (1|Year) + Distance.to.30m, data = TurtleData, family = poisson)

#Initial AIC Scale Selection
summary(NullModel)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: FCAbundance ~ (1 | Year)
##    Data: TurtleData
## 
##      AIC      BIC   logLik deviance df.resid 
##   4975.6   4985.6  -2485.8   4971.6     1078 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4331 -1.0344 -0.2850  0.7805  6.3098 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Year   (Intercept) 0.2325   0.4822  
## Number of obs: 1080, groups:  Year, 10
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.1580     0.1536   7.541 4.67e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
depth.mean.aic <- list("DepthMean250" = DepthMean250,  "DepthMean500" = DepthMean500,  
                       "DepthMean750" = DepthMean750, "DepthMean1000" = DepthMean1000, "Null" = NullModel)

model.sel(depth.mean.aic)
## Model selection table 
##               (Int) MEA.Dpt.250 MEA.Dpt.500 MEA.Dpt.750 MEA.Dpt.1000
## DepthMean1000 1.151                                          -0.1162
## DepthMean750  1.154                            -0.09448             
## DepthMean500  1.156                -0.06834                         
## Null          1.158                                                 
## DepthMean250  1.158    -0.01142                                     
##                     family df    logLik   AICc delta weight
## DepthMean1000 poisson(log)  3 -2459.981 4926.0  0.00      1
## DepthMean750  poisson(log)  3 -2468.583 4943.2 17.20      0
## DepthMean500  poisson(log)  3 -2476.781 4959.6 33.60      0
## Null          poisson(log)  2 -2485.816 4975.6 49.66      0
## DepthMean250  poisson(log)  3 -2485.567 4977.2 51.17      0
## Models ranked by AICc(x) 
## Random terms (all models): 
## '1 | Year'
summary(DepthMean1000)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: FCAbundance ~ (1 | Year) + MEAN.Depth.1000
##    Data: TurtleData
## 
##      AIC      BIC   logLik deviance df.resid 
##   4926.0   4940.9  -2460.0   4920.0     1077 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4152 -1.0250 -0.2146  0.6770  6.5444 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Year   (Intercept) 0.2325   0.4822  
## Number of obs: 1080, groups:  Year, 10
## 
## Fixed effects:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      1.15124    0.15357   7.497 6.55e-14 ***
## MEAN.Depth.1000 -0.11616    0.01617  -7.182 6.86e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## MEAN.D.1000 0.012
slope.mean.aic <- list("SlopeMean250" = SlopeMean250, "SlopeMean500" = SlopeMean500, 
                       "SlopeMean750" = SlopeMean750, "SlopeMean1000" = SlopeMean1000, "Null" = NullModel)

model.sel(slope.mean.aic)
## Model selection table 
##               (Int) MEA.Slp.250 MEA.Slp.500 MEA.Slp.750 MEA.Slp.1000
## SlopeMean250  1.141     -0.1845                                     
## SlopeMean1000 1.148                                          -0.1527
## SlopeMean500  1.149                 -0.1328                         
## SlopeMean750  1.149                              -0.144             
## Null          1.158                                                 
##                     family df    logLik   AICc  delta weight
## SlopeMean250  poisson(log)  3 -2421.862 4849.7   0.00      1
## SlopeMean1000 poisson(log)  3 -2453.637 4913.3  63.55      0
## SlopeMean500  poisson(log)  3 -2454.050 4914.1  64.38      0
## SlopeMean750  poisson(log)  3 -2454.537 4915.1  65.35      0
## Null          poisson(log)  2 -2485.816 4975.6 125.90      0
## Models ranked by AICc(x) 
## Random terms (all models): 
## '1 | Year'
summary(SlopeMean250)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: FCAbundance ~ (1 | Year) + MEAN.Slope.250
##    Data: TurtleData
## 
##      AIC      BIC   logLik deviance df.resid 
##   4849.7   4864.7  -2421.9   4843.7     1077 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3620 -1.0238 -0.1832  0.7159  5.5979 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Year   (Intercept) 0.2325   0.4822  
## Number of obs: 1080, groups:  Year, 10
## 
## Fixed effects:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     1.14113    0.15358    7.43 1.09e-13 ***
## MEAN.Slope.250 -0.18451    0.01642  -11.24  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## MEAN.Sl.250 0.019
rugosity.aic <- list("Rugosity250" = Rugosity250, "Rugosity500" = Rugosity500, 
                      "Rugosity750" = Rugosity750, "Rugosity1000" = Rugosity1000, "Null" = NullModel)

model.sel(rugosity.aic)
## Model selection table 
##              (Int) STD.Dpt.250 STD.Slp.250 STD.Dpt.250:STD.Slp.250 STD.Dpt.500
## Rugosity500  1.158                                                      0.1239
## Rugosity1000 1.181                                                            
## Rugosity250  1.060     0.04737      -0.248                  0.1006            
## Rugosity750  1.168                                                            
## Null         1.158                                                            
##              STD.Slp.500 STD.Dpt.500:STD.Slp.500 STD.Dpt.750 STD.Slp.750
## Rugosity500      -0.3455                -0.05202                        
## Rugosity1000                                                            
## Rugosity250                                                             
## Rugosity750                                          -0.1722      0.1544
## Null                                                                    
##              STD.Dpt.750:STD.Slp.750 STD.Dpt.1000 STD.Slp.1000
## Rugosity500                                                   
## Rugosity1000                              -0.2266       0.2719
## Rugosity250                                                   
## Rugosity750                 -0.05858                          
## Null                                                          
##              STD.Dpt.1000:STD.Slp.1000       family df    logLik   AICc  delta
## Rugosity500                            poisson(log)  5 -2380.883 4771.8   0.00
## Rugosity1000                  -0.08218 poisson(log)  5 -2417.702 4845.5  73.64
## Rugosity250                            poisson(log)  5 -2430.442 4870.9  99.12
## Rugosity750                            poisson(log)  5 -2435.953 4882.0 110.14
## Null                                   poisson(log)  2 -2485.816 4975.6 203.82
##              weight
## Rugosity500       1
## Rugosity1000      0
## Rugosity250       0
## Rugosity750       0
## Null              0
## Models ranked by AICc(x) 
## Random terms (all models): 
## '1 | Year'
summary(Rugosity500)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: FCAbundance ~ (1 | Year) + STD.Depth.500 * STD.Slope.500
##    Data: TurtleData
## 
##      AIC      BIC   logLik deviance df.resid 
##   4771.8   4796.7  -2380.9   4761.8     1075 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5311 -0.9536 -0.2128  0.6746  5.6832 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Year   (Intercept) 0.2325   0.4822  
## Number of obs: 1080, groups:  Year, 10
## 
## Fixed effects:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  1.15750    0.15421   7.506 6.10e-14 ***
## STD.Depth.500                0.12392    0.02882   4.299 1.71e-05 ***
## STD.Slope.500               -0.34553    0.04262  -8.106 5.22e-16 ***
## STD.Depth.500:STD.Slope.500 -0.05202    0.01983  -2.623   0.0087 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) STD.Dp.500 STD.S.
## STD.Dpt.500 -0.027                  
## STD.Slp.500  0.070 -0.728           
## STD.D.500:S -0.086  0.152     -0.511
dist.to.dep.aic <- list("Dist.to.0" = Dist.to.0,
                        "Dist.to.10" = Dist.to.10,
                        "Dist.to.20" = Dist.to.20,
                        "Dist.to.30" = Dist.to.30)

model.sel(dist.to.dep.aic)
## Model selection table 
##            (Int) Dst.to.0m Dst.to.10m Dst.to.20m Dst.to.30m       family df
## Dist.to.30 1.123                                     0.2686 poisson(log)  3
## Dist.to.0  1.125   -0.2583                                  poisson(log)  3
## Dist.to.20 1.126                          0.2517            poisson(log)  3
## Dist.to.10 1.158              0.03017                       poisson(log)  3
##               logLik   AICc  delta weight
## Dist.to.30 -2357.706 4721.4   0.00  0.917
## Dist.to.0  -2360.435 4726.9   5.46  0.060
## Dist.to.20 -2361.376 4728.8   7.34  0.023
## Dist.to.10 -2484.016 4974.1 252.62  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
## '1 | Year'
summary(Dist.to.20)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: FCAbundance ~ (1 | Year) + Distance.to.20m
##    Data: TurtleData
## 
##      AIC      BIC   logLik deviance df.resid 
##   4728.8   4743.7  -2361.4   4722.8     1077 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3598 -0.9920 -0.1704  0.7198  5.5644 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Year   (Intercept) 0.2325   0.4822  
## Number of obs: 1080, groups:  Year, 10
## 
## Fixed effects:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      1.12581    0.15361   7.329 2.32e-13 ***
## Distance.to.20m  0.25169    0.01594  15.790  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## Dstnc.t.20m -0.027
Turt.FC.data1 <- data.frame(FCs = TurtleData$FCAbundance,
                         Year = as.numeric(TurtleData$Year),
                         muDep = abs(TurtleData$MEAN.Depth.1000),
                         muSlop = TurtleData$MEAN.Slope.250,
                         Rugosity = TurtleData$STD.Slope.500,
                         dist30 = TurtleData$Distance.to.30m,
                         shelfW = (TurtleData$Distance.to.20m - TurtleData$Distance.to.10m)
                         )



ggpairs(Turt.FC.data1, columns = 3:7) + theme_bw()

GlobalModel <- glmer(FCs ~ muDep + muSlop + Rugosity + 
             dist30 + shelfW + (1|Year), data = Turt.FC.data1, family = poisson)

performance::check_collinearity(GlobalModel)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##      Term  VIF Increased SE Tolerance
##     muDep 1.08         1.04      0.93
##    muSlop 1.74         1.32      0.58
##  Rugosity 1.60         1.26      0.63
##    dist30 2.66         1.63      0.38
##    shelfW 1.52         1.23      0.66
GlobalModel <- glmer(FCs ~ muDep + muSlop + Rugosity + 
             shelfW + (1|Year), data = Turt.FC.data1, family = poisson)

performance::check_collinearity(GlobalModel)
## # Check for Multicollinearity
## 
## Low Correlation
## 
##      Term  VIF Increased SE Tolerance
##     muDep 1.07         1.04      0.93
##    muSlop 1.09         1.04      0.92
##  Rugosity 1.31         1.14      0.76
##    shelfW 1.25         1.12      0.80
dredge(GlobalModel)
## Fixed term is "(Intercept)"
## Global model call: glmer(formula = FCs ~ muDep + muSlop + Rugosity + shelfW + (1 | 
##     Year), data = Turt.FC.data1, family = poisson)
## ---
## Model selection table 
##    (Intrc)    muDep   muSlp   Rgsty  shlfW df    logLik   AICc  delta weight
## 16   1.168 -0.06482 -0.1502 -0.1566 0.1285  6 -2323.531 4659.1   0.00  0.657
## 15   1.113          -0.1542 -0.1475 0.1352  5 -2325.194 4660.4   1.30  0.343
## 11   1.121          -0.1804         0.1817  4 -2344.831 4697.7  38.56  0.000
## 12   1.147 -0.02961 -0.1793         0.1796  5 -2344.470 4699.0  39.85  0.000
## 8    1.237 -0.13510 -0.1311 -0.2352         5 -2353.371 4716.8  57.66  0.000
## 7    1.122          -0.1373 -0.2233         4 -2361.025 4730.1  70.95  0.000
## 14   1.217 -0.11230         -0.2147 0.1155  5 -2361.194 4732.4  73.30  0.000
## 13   1.122                  -0.1998 0.1272  4 -2366.224 4740.5  81.34  0.000
## 6    1.272 -0.16910         -0.2843         4 -2383.573 4775.2 116.04  0.000
## 5    1.129                  -0.2707         3 -2395.585 4797.2 138.05  0.000
## 10   1.197 -0.07077                 0.1929  4 -2399.262 4806.6 147.42  0.000
## 9    1.135                          0.1979  3 -2401.327 4808.7 149.54  0.000
## 4    1.245 -0.12060 -0.1809                 4 -2415.528 4839.1 179.95  0.000
## 3    1.141          -0.1845                 3 -2421.862 4849.7 190.60  0.000
## 2    1.296 -0.15990                         3 -2474.804 4955.6 296.49  0.000
## 1    1.158                                  2 -2485.816 4975.6 316.50  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
## '1 | Year'
res <- dredge(GlobalModel)
## Fixed term is "(Intercept)"
importance(res)
##                      muSlop shelfW Rugosity muDep
## Sum of weights:      1.00   1.00   1.00     0.66 
## N containing models:    8      8      8        8
car::Anova(GlobalModel)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: FCs
##           Chisq Df Pr(>Chisq)    
## muDep     3.327  1    0.06815 .  
## muSlop   73.968  1  < 2.2e-16 ***
## Rugosity 38.080  1  6.791e-10 ***
## shelfW   58.317  1  2.231e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(Turt.FC.data1, aes(muDep, FCs)) +
  geom_point(size = 1) +
   geom_smooth(method = "glm", se = FALSE,
       method.args = list(family = "poisson"))
## `geom_smooth()` using formula 'y ~ x'

ggplot(Turt.FC.data1, aes(muSlop, FCs)) +
  geom_point(size = 1) +
   geom_smooth(method = "glm", se = FALSE,
       method.args = list(family = "poisson"))
## `geom_smooth()` using formula 'y ~ x'

ggplot(Turt.FC.data1, aes(Rugosity, FCs)) +
  geom_point(size = 1) +
   geom_smooth(method = "glm", se = FALSE,
       method.args = list(family = "poisson"))
## `geom_smooth()` using formula 'y ~ x'

ggplot(Turt.FC.data1, aes(shelfW, FCs)) +
  geom_point(size = 1) +
   geom_smooth(method = "glm", se = FALSE,
       method.args = list(family = "poisson"))
## `geom_smooth()` using formula 'y ~ x'